
!-------------------------------------------------------------------------------------
!
!    call evolve_F_long_advect( dfu, dt)
!
!    call evolve_F_trans_advect( dfu, Fbgk, uxMat, dt)
!
!    call evolve_efield( dfu, Ex2D, Ez2D)
!
!    call evolve_F_accel( dfu, Ebgk, DVbgk, Ez2D, dt)
!
!
!   call evolve_U_long_advect( uxMat, dt)
!
!   call evolve_U_trans_advect( uxMat, dt)
!
!   call evolve_U_sourceE( uxMat, Ex2D, dt)



!-------------------------------------------------------------------------------------
subroutine evolve_lin_F_long_advect(dfu, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, Dzzzz, SideLoss
  use fft_wrap
  use fft_3d

  implicit none

  complex, intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: dt

  integer                :: iz, ix, iy, iDir
  complex                :: fz(jz)
  complex                :: iVdt(jv), exp_iVdt(jv)

  !-- advect in z-direction

  call dfftw_execute_dft(plan_z5forward, dfu, dfu)
  dfu = dfu/jz

  fz(:) = Dzzzz * kZvec(:)**4 *dt + SideLoss*dt
  iVdt(:) = iunit * dt * vVec(:)

  !OMP PARALLEL PRIVATE (iv,iVdt) 
  !OMP DO 


  do iDir=1,jDir
     do ix=1, jx
        do iy=1, jy

           !$OMP PARALLEL PRIVATE (iz,exp_iVdt) 
           !$OMP DO 
           do iz=1,jz

              exp_iVdt(:) =  exp( - iVdt(:) * kZvec(iz) - fz(iz) )

              dfu(:, iz, iy, ix, iDir) = dfu(:, iz, iy, ix, iDir) * exp_iVdt(:)

           enddo
           !$OMP END DO
           !$OMP END PARALLEL

        enddo
     enddo
  enddo


  !OMP END DO
  !OMP END PARALLEL

  call dfftw_execute_dft(plan_z5backward, dfu, dfu)

end subroutine evolve_lin_F_long_advect


!-------------------------------------------------------------------------------------
subroutine evolve_lin_F_trans_advect(dfu, Fbgk, uxMat, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, ux0
  use fft_wrap
  use fft_3d

  implicit none

  complex, intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: Fbgk(jv, jz)
  complex, intent(inout) :: uxMat(jz, jy, jx, jDir)   ! <-- should be just intent(in)
  real,    intent(in)    :: dt

  complex                :: c, e
  integer                :: ix, iz, iy, iv, iDir


  !-- advect in x-direction

  call dfftw_execute_dft(plan_x5forward, dfu, dfu)
  dfu = dfu/jx

  call dfftw_execute_dft(plan_x4forward, uxMat, uxMat)
  uxMat = uxMat/jx

  do iDir=1,jDir
     do ix=1,jx
        do iy=1,jy

           e = -iunit * ux0(iDir) *  kXvec(ix)
           e =  exp(e*dt) 

           !$OMP PARALLEL PRIVATE(iz, c)
           !$OMP DO
           do iz=1,jz
              c = uxMat(iz,iy,ix,iDir) / ux0(iDir)
              dfu(:, iz, iy, ix, iDir) =  dfu(:, iz, iy, ix, iDir)*e  +  c*Fbgk(:,iz)*(e - 1)
           enddo
           !$OMP END DO
           !$OMP END PARALLEL

        enddo
     enddo
  enddo

  call dfftw_execute_dft(plan_x4backward, uxMat, uxMat)

  call dfftw_execute_dft(plan_x5backward, dfu, dfu)


end subroutine evolve_lin_F_trans_advect


!-------------------------------------------------------------------------------------
subroutine evolve_lin_F_accel(dfu, Ebgk, ftDVbgk, Ez2D, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, Dvvvv, CollisionRate
  use fft_wrap
  use fft_3d

  implicit none

  complex, intent(inout) :: dfu(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: Ebgk(jz)
  complex, intent(in)    :: ftDVbgk(jv, jz)
  complex, intent(in)    :: Ez2D(jz, jy, jx)
  real,    intent(in)    :: dt

  integer                :: iz, iv, ix, iy, iDir
  complex                :: b(jv), c(jv), e(jv), fv(jv), DvvvvTerm(jv), b0(jv)



  DvvvvTerm(:) =  Dvvvv*kVvec(:)**4

  b0(:)  =  DvvvvTerm(:) + CollisionRate * kVvec(:)**2
  
  !-- advect in v-direction

  call dfftw_execute_dft(plan_v5forward, dfu, dfu)
  dfu = dfu/jv

   do idir=1,jDir
     do ix=1, jx
     do iy=1, jy

        !-- before applying hyperviscosity, subtract z-average, on Denis's demand
        fv(:) = sum(dfu(:,:,iy,ix,iDir), 2)/jz


        !$OMP PARALLEL PRIVATE(iz, b, e, c)
        !$OMP DO

        do iz=1,jz

           b(:) = - iunit*Ebgk(iz)*kVvec(:) - b0(:)
           e(:) = exp(b*dt)

           c(:) = - ftDVbgk(:,iz)*Ez2D(iz,iy,ix)  + DvvvvTerm * fv
           c(:) = c/b
           c(1) = 0

           dfu(:, iz, iy, ix, iDir) = dfu(:, iz, iy, ix, iDir) * e + c*(e-1)
    
        end do

        !$OMP END DO
        !$OMP END PARALLEL

     end do
     end do
  end do

  call dfftw_execute_dft(plan_v5backward, dfu, dfu)

 
end subroutine evolve_lin_F_accel


!*************************************************************************************

!-------------------------------------------------------------------------------------
subroutine rk_rhs(s, Ubgk, Ex2D)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, vPhase, ux0
  use fft_wrap
  use fft_3d

  implicit none

  complex, intent(inout) :: s(jz, jy, jx, jdir)
  complex, intent(inout) :: Ex2D(jz, jy, jx)
  real,    intent(in)    :: Ubgk(jz)
 
  integer                :: iv, ix, iy, iDir

  complex                :: s0(jz, jy, jx, jdir)

  s0=s;

  call dfftw_execute_dft(plan_z4forward, s, s)
  s = s/jz

  do iDir=1,jDir
     do ix=1,jx
        do iy=1,jy
            s(:,iy,ix,iDir) = - s(:,iy,ix,iDir) * iunit * kZvec(:) 
       enddo
     enddo
  enddo
  
  call dfftw_execute_dft(plan_z4backward, s, s)

  do iDir=1,jDir
     do ix=1,jx
        do iy=1,jy          
           !s(:,iy,ix,iDir) = s(:,iy,ix,iDir)*Ubgk(:) - iunit*kXvec(1)*ux0(iDir) * s0(:,iy,ix,iDir) + Ex2D(:,iy,ix)
           s(:,iy,ix,iDir) = s(:,iy,ix,iDir)*(-vPhase) - iunit*kXvec(1)*ux0(iDir) * s0(:,iy,ix,iDir)  + Ex2D(:,iy,ix)
        end do
     end do
  end do

end subroutine rk_rhs


!-------------------------------------------------------------------------------------
subroutine evolve_lin_U_rk(uxMat, Ex2D, Ubgk, dt)
!-------------------------------------------------------------------------------------

  use dimensions

  implicit none

  complex, intent(inout) :: uxMat(jz, jy, jx, jdir)
  complex, intent(inout) :: Ex2D(jz, jy, jx)
  real,    intent(in)    :: Ubgk(jz)
  real,    intent(in)    :: dt

  complex                :: s(jz, jy, jx, jdir)
  complex                :: s0(jz, jy, jx, jdir)


  s = uxMat;
  call rk_rhs(s, Ubgk, Ex2D)
  s0 = s

  s = uxMat + s*dt/2
  call rk_rhs(s, Ubgk, Ex2D)
  s0 = s0 + 2*s

  s = uxMat + s*dt/2
  call rk_rhs(s, Ubgk, Ex2D)
  s0 = s0 + 2*s

  s = uxMat + s*dt
  call rk_rhs(s, Ubgk, Ex2D)
  s0 = s0 + s

  uxMat = uxMat + s0*dt/6

end subroutine evolve_lin_U_rk

!-------------------------------------------------------------------------------------
subroutine evolve_lin_U (uxMat, Ex2D, Ubgk, dt)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only: iunit, vPhase, iWaveFrame, ux0
  use fft_3d
  use fft_wrap

  implicit none

  complex, intent(inout) :: uxMat(jz, jy, jx, jdir)
  complex, intent(inout) :: Ex2D(jz, jy, jx)
  real,    intent(in)    :: Ubgk(jz)
  real,    intent(in)    :: dt

  integer                :: iv, ix, iy, iDir
  complex                :: c(jz), b(jz), e(jz)


  call dfftw_execute_dft(plan_z4forward, uxMat, uxMat)
  uxMat = uxMat/jz

  call dfftw_execute_dft(plan_xyzforward, Ex2D, Ex2D)
  Ex2D = Ex2D/(jx*jy*jz)

  do iDir=1,jDir
     do iy=1,jy
        do ix=1,jx

           b = - iunit*ux0(iDir)*kXvec(ix)  -  iunit*(-vPhase*iWaveFrame)*kZvec(:)
           e = exp(b*dt)
           c = Ex2D(:,iy,ix)/b
           if (b(1) .eq. 0) c(1)=0

           !$OMP PARALLEL SHARED(uxMat,c,e)
           !$OMP WORKSHARE
           uxMat(:,iy,ix,iDir) = uxMat(:,iy,ix,iDir) * e + c*(e-1)
           !$OMP END WORKSHARE NOWAIT
           !$OMP END PARALLEL

        enddo
     enddo
  enddo


  call dfftw_execute_dft(plan_z4backward, uxMat, uxMat)

  call dfftw_execute_dft(plan_xyzbackward, Ex2D, Ex2D)


end subroutine evolve_lin_U



!-------------------------------------------------------------------------------------
subroutine mult_4D(u, m)
!-------------------------------------------------------------------------------------

  use dimensions

  implicit none

  complex, intent(inout) :: u(jz, jy, jx, jDir)
  real,    intent(in)    :: m

  !$OMP PARALLEL SHARED(u)
  !$OMP WORKSHARE
  u = u * m
  !$OMP END WORKSHARE NOWAIT
  !$OMP END PARALLEL
  
end subroutine mult_4D

!-------------------------------------------------------------------------------------
subroutine mult_5D(u, m)
!-------------------------------------------------------------------------------------

  use dimensions

  implicit none

  complex, intent(inout) :: u(jv, jz, jy, jx, jDir)
  real,    intent(in)    :: m

  !$OMP PARALLEL SHARED(u)
  !$OMP WORKSHARE
  u = u * m
  !$OMP END WORKSHARE NOWAIT
  !$OMP END PARALLEL
  
end subroutine mult_5D



!-------------------------------------------------------------------------------------
